library(tidyverse)
library(ggpubr)
library(scales)
library(cowplot)


## read in data for figure 1

brfss_1 <- read_csv("dis_effects_brfss.csv")
hps_1 <- read_csv("dis_effects_hps.csv")
ncvs_1 <- read_csv("ncvs dis_effects.csv")

# rename hps variables for consistency with brfss and ncvs
hps_1$phys_dis_fac <- hps_1$dis_phys_fac

## make figure 1

p1 <- brfss_1 %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(plot.title = element_text(hjust = 0.5))

p2 <- hps_1 %>%
  filter(TBIRTH_YEAR >= 1940) %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5))

p3 <- ncvs_1 %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "NCVS",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5))

annotate_figure(ggarrange(p1, p2, p3, ncol = 3, legend = "bottom", common.legend = T))

## Figure 2

diabetes <- read_csv("diabetes_effects.csv")

p10 <- diabetes %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = diabetes_class, linetype = diabetes_class, fill = diabetes_class)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Diabetes type",
       color = "Diabetes type",
       linetype = "Diabetes type") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5))

## figure 3 is created in Combined Analysis.R

## read in data for figure 4

brfss_afab <- read_csv("afab_effects.csv")
hps_afab <- read_csv("dis_effects_afab_hps.csv")
ncvs_afab <- read_csv("ncvs dis_effects_afab.csv")

brfss_amab <- read_csv("amab_effects_brfss.csv")
hps_amab <- read_csv("dis_effects_amab_hps.csv")
ncvs_amab <- read_csv("ncvs dis_effects_amab.csv")

# rename hps variables for consistency with brfss and ncvs
hps_afab$phys_dis_fac <- hps_afab$dis_phys_fac
hps_amab$phys_dis_fac <- hps_amab$dis_phys_fac


p4 <- brfss_afab %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned female at birth (BRFSS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

p5 <- hps_afab %>%
  filter(TBIRTH_YEAR >= 1940) %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned female at birth (HPS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

p6 <- ncvs_afab %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned female at birth (NCVS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

p7 <- brfss_amab %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned male at birth (BRFSS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

p8 <- hps_amab %>%
  filter(TBIRTH_YEAR >= 1940) %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned male at birth (HPS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

p9 <- ncvs_amab %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = phys_dis_fac, linetype = phys_dis_fac, fill = phys_dis_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned male at birth (NCVS)",
       fill = "Physical disability",
       color = "Physical disability",
       linetype = "Physical disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.15), oob = rescale_none) +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

annotate_figure(ggarrange(p4, p5, p6, p7, p8, p9, nrow = 2, ncol = 3, common.legend = T, legend = "bottom"))

## Figure 5

emmeans_phys_df <- read_csv("emmeans_phys_df.csv")

plot1 <- emmeans_phys_df %>% # physical disability
  ggplot(aes(x = prop_dis, y = prob, col = phys_dis, linetype = phys_dis, fill = phys_dis, ymin = lower, ymax = upper)) +
  geom_line() +
  geom_ribbon(alpha = 0.3) +
  geom_rug(data = regression_vars, mapping = aes(x = prop_dis, y = NULL, col = NULL,
                                                 ymax = NULL, ymin = NULL), sides = "b") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.25), oob = rescale_none) +
  scale_x_continuous(limits = c(0.02, 0.25)) +
  labs(x = "Proportion of Students with Physical Disabilities",
       y = "Probability of Identifying as Trans",
       linetype = "Physical Disability Status",
       col = "Physical Disability Status",
       fill = "Physical Disability Status",
       title = "Gender Identification by School Physical Disability Composition") +
  scale_color_hue(labels = c("No", "Yes")) +
  scale_fill_hue(labels = c("No", "Yes")) +
  scale_linetype_discrete(labels = c("No", "Yes"))  +
  theme(legend.position = "bottom", , plot.title = element_text(hjust = 0.5))

plot2 <- emmeans_ment_df %>% # mental disability
  ggplot(aes(x = prop_mdis, y = prob, col = ment_dis_fac, linetype = ment_dis_fac, fill = ment_dis_fac, ymin = lower, ymax = upper)) +
  geom_line() +
  geom_ribbon(alpha = 0.3) +
  geom_rug(data = regression_vars, mapping = aes(x = prop_mdis, y = NULL, col = NULL,
                                                 ymax = NULL, ymin = NULL), sides = "b") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.25), oob = rescale_none) +
  scale_x_continuous(limits = c(0.05, 0.4)) +
  labs(x = "Proportion of Students with Mental Disabilities",
       y = "Probability of Identifying as Trans",
       linetype = "Mental Disability Status",
       col = "Mental Disability Status",
       fill = "Mental Disability Status",
       title = "Gender Identification by School Mental Disability Composition") +
  scale_color_hue(labels = c("No", "Yes")) +
  scale_fill_hue(labels = c("No", "Yes")) +
  scale_linetype_discrete(labels = c("No", "Yes"))  +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

plot1 + plot2 # Figure 5


## Tables 2-4
brfss_wide <- pivot_wider(brfss_1, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

brfss_wide <- brfss_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

brfss_wide

hps_wide <- pivot_wider(hps_1, id_cols = c("TBIRTH_YEAR"), names_from = "dis_phys_fac", values_from = "est")

hps_wide <- hps_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

hps_wide


ncvs_wide <- pivot_wider(ncvs_1, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

ncvs_wide <- ncvs_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

ncvs_wide

brfss_afab_wide <- pivot_wider(brfss_afab, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

brfss_afab_wide <- brfss_afab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

brfss_afab_wide

hps_afab_wide <- pivot_wider(hps_afab, id_cols = c("TBIRTH_YEAR"), names_from = "dis_phys_fac", values_from = "est")

hps_afab_wide <- hps_afab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

hps_afab_wide

ncvs_afab_wide <- pivot_wider(ncvs_afab, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

ncvs_afab_wide <- ncvs_afab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)
ncvs_afab_wide

brfss_amab_wide <- pivot_wider(brfss_amab, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

brfss_amab_wide <- brfss_amab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)
         

brfss_amab_wide

hps_amab_wide <- pivot_wider(hps_amab, id_cols = c("TBIRTH_YEAR"), names_from = "dis_phys_fac", values_from = "est")

hps_amab_wide <- hps_amab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

hps_amab_wide


ncvs_amab_wide <- pivot_wider(ncvs_amab, id_cols = c("birth_year"), names_from = "phys_dis_fac", values_from = "est")

ncvs_amab_wide <- ncvs_amab_wide %>%
  mutate(risk_ratio = Yes/No,
         marginal_effect = Yes - No)

ncvs_amab_wide

## Appendix B

app_b_all <- read_csv("dis_effects_hps_app_b.csv")
app_b_afab <- read_csv("dis_effects_afab_hps_app_b.csv")
app_b_amab <- read_csv("dis_effects_amab_hps_app_b.csv")

b1_1 <- app_b_all %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = dis_phys_fac, linetype = dis_phys_fac, fill = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       fill = "Physical Disability",
       color = "Physical Disability",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom")

b1_afab <- app_b_afab %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = dis_phys_fac, linetype = dis_phys_fac, fill = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned Female at Birth",
       fill = "Physical Disability",
       color = "Physical Disability",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

b1_amab <- app_b_amab %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = dis_phys_fac, linetype = dis_phys_fac, fill = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned Male at Birth",
       fill = "Physical Disability",
       color = "Physical Disability",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

library(ggpubr)

b1_1
ggarrange(b1_afab, b1_amab, common.legend = T, legend = "bottom")

##appendix B black and white

b1_bw <- app_b_all %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom")

b1_afab_bw <- app_b_afab %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned Female at Birth",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

b1_amab_bw <- app_b_amab %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = dis_phys_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Assigned Male at Birth",
       linetype = "Physical Disability") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.20), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

b1_bw
ggarrange(b1_afab_bw, b1_amab_bw, common.legend = T, legend = "bottom")


## Appendix C

brfss_blind <- read_csv("brfss_blind.csv")
brfss_deaf <- read_csv("brfss_deaf.csv")
brfss_walk <- read_csv("brfss_walk.csv")
brfss_dress <- read_csv("brfss_dress.csv")

hps_blind <- read_csv("hps_blind_effects.csv")
hps_deaf <- read_csv("hps_deaf_effects.csv")
hps_walk <- read_csv("hps_walk_effects.csv")
hps_dress <- read_csv("hps_dress_effects.csv")

c1_1 <- brfss_blind %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = blind_fac, linetype = blind_fac, fill = blind_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Blind/Difficulty Seeing",
       color = "Blind/Difficulty Seeing",
       linetype = "Blind/Difficulty Seeing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_2 <- hps_blind %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = blind_fac, linetype = blind_fac, fill = blind_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Blind/Difficulty Seeing",
       color = "Blind/Difficulty Seeing",
       linetype = "Blind/Difficulty Seeing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_3 <- brfss_deaf %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = deaf_fac, linetype = deaf_fac, fill = deaf_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Deaf/Difficulty Hearing",
       color = "Deaf/Difficulty Hearing",
       linetype = "Deaf/Difficulty Hearing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_4 <- hps_deaf %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = deaf_fac, linetype = deaf_fac, fill = deaf_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Deaf/Difficulty Hearing",
       color = "Deaf/Difficulty Hearing",
       linetype = "Deaf/Difficulty Hearing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_5 <- brfss_walk %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = walk_fac, linetype = walk_fac, fill = walk_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Difficulty Walking/Climbing Stairs",
       color = "Difficulty Walking/Climbing Stairs",
       linetype = "Difficulty Walking/Climbing Stairs") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_6 <- hps_walk %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = walk_fac, linetype = walk_fac, fill = walk_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Difficulty Walking/Climbing Stairs",
       color = "Difficulty Walking/Climbing Stairs",
       linetype = "Difficulty Walking/Climbing Stairs") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_7 <- brfss_dress %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, col = dress_fac, linetype = dress_fac, fill = dress_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       fill = "Difficulty Dressing/Bathing",
       color = "Difficulty Dressing/Bathing",
       linetype = "Difficulty Dressing/Bathing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_8 <- hps_dress %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = dress_fac, linetype = dress_fac, fill = dress_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Difficulty Dressing/Bathing",
       color = "Difficulty Dressing/Bathing",
       linetype = "Difficulty Dressing/Bathing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

ggarrange(c1_1, c1_3, c1_5, c1_7,
          c1_2, c1_4, c1_6, c1_8,
          ncol = 4, nrow = 2)

## Appendix C black and white
c1_1_bw <- brfss_blind %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, linetype = blind_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       linetype = "Blind/Difficulty Seeing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_2_bw <- hps_blind %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = blind_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       fill = "Blind/Difficulty Seeing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_3_bw <- brfss_deaf %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, linetype = deaf_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       linetype = "Deaf/Difficulty Hearing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_4_bw <- hps_deaf %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = deaf_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       linetype = "Deaf/Difficulty Hearing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_5_bw <- brfss_walk %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, linetype = walk_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       linetype = "Difficulty Walking/Climbing Stairs") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_6_bw <- hps_walk %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = walk_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       linetype = "Difficulty Walking/Climbing Stairs") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_7_bw <- brfss_dress %>%
  filter(birth_year >= 1940) %>%
  ggplot(aes(x = birth_year, y = est, linetype = dress_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "BRFSS",
       linetype = "Difficulty Dressing/Bathing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

c1_8_bw <- hps_dress %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = dress_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "HPS",
       linetype = "Difficulty Dressing/Bathing") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.33), oob = rescale_none) +
  theme(legend.position = "bottom")

ggarrange(c1_1_bw, c1_3_bw, c1_5_bw, c1_7_bw,
          c1_2_bw, c1_4_bw, c1_6_bw, c1_8_bw,
          ncol = 4, nrow = 2)

## Appendix D
hps_severe <- read_csv("dis_effects_severe.csv")

hps_severe$dis_severe_fac <- factor(hps_severe$dis_severe_fac,
                                    levels = c("No difficulty", "Some difficulty", "A lot of difficulty",
                                               "Cannot do at all"),
                                    ordered = T)

d1 <- hps_severe %>%
  filter(TBIRTH_YEAR > 1935) %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, col = dis_severe_fac, linetype = dis_severe_fac, fill = dis_severe_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Disability Severity (HPS)",
       fill = "Degree of Impairment",
       color = "Degree of Impairment",
       linetype = "Degree of Impairment") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.5), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))


d1_bw <- hps_severe %>%
  filter(TBIRTH_YEAR > 1935) %>%
  ggplot(aes(x = TBIRTH_YEAR, y = est, linetype = dis_severe_fac)) +
  geom_ribbon(aes(ymin = est - 1.96*se,
                  ymax = est + 1.96*se), alpha = 0.3) +
  geom_path() +
  labs(x = "Birth year",
       y = "Probability of identifying as transgender",
       title = "Disability Severity (HPS)",
       linetype = "Degree of Impairment") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, 0.5), oob = rescale_none) +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
